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ABSTRACT 

The origin of the Fermi bubbles recently detected by the Fermi Gamma-ray Space Telescope in the 
inner Galaxy is mysterious. In the companion paper Guo & Mathews (Paper I), we use hydrodynamic 
simulations to show that they could be produced by a recent powerful AGN jet event. Here we further 
explore this scenario to study the potential roles of shear viscosity and cosmic ray (CR) diffusion on the 
morphology and CR distribution of the bubbles. We show that even a relatively low level of viscosity 
(^visc ^ 3 g cm~^ s~^, or ~ 0.1% - 1% of Braginskii viscosity in this context) could effectively suppress 
the development of Kelvin-Helmholtz instabilities at the bubble surface, resulting in smooth bubble 
edges as observed. Furthermore, viscosity reduces circulating motions within the bubbles, which would 
otherwise mix the CR-carrying jet backflow near bubble edges with the bubble interior. Thus viscosity 
naturally produces an edge-favored CR distribution, an important ingredient to produce the observed 
flat gamma-ray surface brightness distribution. Generically, such a CR distribution often produces 
a limb-brightened gamma-ray intensity distribution. However, we show that by incorporating CR 
diffusion which is strongly suppressed across the bubble surface (as inferred from sharp bubble edges) 
but is close to canonical values in the bubble interior, we obtain a reasonably flat gamma-ray intensity 
profile. The similarity of the resulting CR bubble with the observed Fermi bubbles strengthens our 
previous result in Paper I that the Fermi bubbles were produced by a recent AGN jet event. Studies 
of the nearby Fermi bubbles may provide a unique opportunity to study the potential roles of plasma 
viscosity and CR diffusion on the evolution of AGN jets and bubbles. 

Subject headings: cosmic rays - galaxies: active - galaxies: jets - Galaxy: nucleus - gamma rays: 
galaxies 
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1. INTRODUCTION 

The recent discovery of two large gamm a-ray bubbles 
by th e Fermi Gamma -Ray Space Telescope (jDobler et al.l 
I2Q1QI : iSu et al.l l2QlQf ) has significantly changed the big 
picture of our Galaxy, the Milky Way. These ''Fermi 
bubbles" emit at 1 < < 100 GeV in the inner Galaxy, 
are nearly symmetric about the Galactic plane, and ex- 
tend to ^ 50° (~ 10 kpc) above and below the Galactic 
center (GC), with a width of about 40° in longitude. In 
addition, they have approximately uniform gamma-ray 
surface brightness with sharp edges. 

The origin of the Fermi bubbles has recently received 
a lot of attention and is actively debated in the lit- 
erature. The sharp edges and bilobular morphology 
of the bubbles make them difficult to be explained by 
either diffused CRs from the Galactic di sk or annihi- 
lations of dark matter particles (though iDobler et al.l 
12011. claimed that dark matter annihilations in a pro- 
late halo with a nisotropic CR diffusion may explain the 
latter feature). iCrocker fc AharonianI (|2011f ) suggested 
that the gamma ray emission is powered by CR pro- 
tons, which are continuously injected by supernova ex- 
plosions in the Galactic center during a few Gyrs, though 
this seems difficult to reconcile with the sharp edges 
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which implies a mo re re cent trans i ent ev ent. In contrast, 
IDobler et all (j2010f ) and ISu et all (|2010^ argued that the 
emission may be dominated by upscattering of photons 
in the interstellar radiation field (ISRF) and the cos- 
mic microwave background by CR electrons, whose syn- 
chrotron emission may have already been detected at tens 
of GHz b y the Wilkinson Microwave Anisotropy Probe 
fWMAP: iFinkbeinerl |2QQ1 Dobler & Finkbeiner 2008|). 
iCheng et al.l (|2011f ) argued that periodic star capture 
processes by the Galactic supermassive black hole, Sgr 
A*, release AGN winds into the Galactic halo, produc- 
i ng th e Fermi bubbles within a few Myrs. IZubovas et al.l 
(|2011[ ) suggested that the near-spherical outflow from a 
quasar event of Sgr A* around 6 Myr ago may explain 

the origin of the bubbles. 

In a companion paper (|Guo fc Mathewsl 120121 : here- 
after denoted as "Paper I"), we performed the first nu- 
merical simulation following the dynamical evolution of 
the Fermi bubbles, and showed that a recent AGN jet 
event originated from Sgr A* around 1-3 Myr ago 
can reproduce the Fermi bubbles with roughly the ob- 
served location, size, and shape. Our jet model is in- 
spired by many extragalactic AGN jets, which are clearly 
produ cing CR-filled bubbles seen in radio observations 
(McN amara fc NulsenI l2QQ7l ) . In our model, the oppos- 
ing jets, dominated by kinetic energy and over-pressured 
by either CR or thermal pressure, were active for ~ 0.1 
- 0.5 Myr and moderately light. We also show that the 
sharp bubble edges require that CR diffusion across the 
bubble edges is suppressed significantly below the CR 
diffusion rate estimated in the solar vicinity. 
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While many observational features of the Fermi bub- 
bles, particularly their age, location, size and shape, are 
reproduced by our model in Paper I, the simulated bub- 
ble is deficient in two important respects - surface ir- 
regularities and limb darkening in gamma ray surface 
brightness, disagreeing significantly with smooth edges 
and roughly uniform gamma ray surface brightness of 
the observed Fermi bubbles. However, these inconsisten- 
cies do not necessarily mean that the jet scenario for the 
Fermi bubbles is wrong. Instead, they may be smoking- 
gun signatures of additional physics, which plays a sig- 
nificant role during the jet evolution. Surface irregu- 
larities induced by Kelvin-Helmholtz instabilities have 
previously been seen in hydrodynamical simulations of 
buoyantly-rising X-ray cavities in galaxy clusters, where 
add itional physical mechanisms, including hot gas viscos- 
ity (|Revnolds et a l. 2005; Kaiser et al. 20 Q5D and mag- 
netic draping (Lvutikov 2006; Ruszkowski et al.l 120071: 
iDursi fc Pfrommerll2QQi ).^ have been invoked to suppress 
the instabilities. The fiatness of the line-of-sight pro- 
jected gamma ray intensity distribution is not trivial, as 
a spatially-uniform CR distribution produces a center- 
brightened gamma ray surface brightness while an edge- 
dominated CR distribution produces limb brightening. 
The difficulties involved in producing such a 'ffat' gamma 
ray brig htness in the Ferm i bubb les were previously no- 
ticed by iMertsch fc Sarkarl ()2011l ) , who suggested that it 
may be achieved if CR electrons are re-accelerated pref- 
erentially near bubble edges. 

In this paper, we further explore our jet model for the 
Fermi bubbles developed in Paper I by including ad- 
ditional gas microphysics - shear viscosity. We show 
that even a relatively low level of viscosity (compared 
to Spitzer viscosity in the shock-heated surrounding gas) 
significantly affects the evolution of the resulting Fermi 
bubbles, helping produce smooth bubble edges and a 
ffat gamma ray surface brightness as observed. We in- 
vestigate the roles of viscosity and CR diffusion on the 
morphology and CR distribution of the bubbles by di- 
rectly follo wing the jet evolution , which differs signiff- 
cantly from iRevnolds et al.l (|2005l ) , who investigated the 
role of viscosity on the buoyant rise of initially static 
bubbles in galaxy clusters. Furthermore, our paper is 
the ffrst to study the role of viscosity in the context of 
the Fermi bubbles. 

The rest of the paper is organized as follows. In Sec- 
tion O we describe our model and numerical methods. 
We present our results in Section [H and summarize our 
main conclusions with implications in Section HI In the 
Appendix, we explicitly present how we implement shear 
viscosity in cylindrical coordinates with axisymmetry. 

2. THE MODEL AND NUMERICAL METHODS 

In Paper I and this paper, we study the formation of 
the Fermi bubbles in the Milky Way's potential well with 
AGN jets using numerical simulations. We presented the 
basic picture of the jet scenario in Paper I, where we 
show that the bubbles can be formed by a recent AGN 
jet event which started around 1-3 Myr ago and lasted 
for ~ 0.1 - 0.5 Myr. In the current paper, we continue 
to investigate the potential roles of gas viscosity and CR 
diffusion on the evolution of the Fermi bubbles in the jet 
scenario, while directly comparing the model to Fermi 
observations. We refer the reader to Paper I for details of 
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Fig. 1. — Top: The initial density distribution of the hot ther- 
mal gas in the Galactic halo for all runs presented in this paper. 
Horizontal and vertical axes refer to R and z respectively, labeled 
in kpc. The hot gas is assumed to be isothermal (T = 2.4 x 10^ 
K) and in hydrostatic equilibrium at time t = 0. Bottom: The 
density distribution of the hot gas at t = tpermi in run V3-diff3 
(see Table 1). The low-density cavity and the surrounding shock, 
both produced by the AGN jet event, are clearly seen. 

our model and assumptions. Here we simply summarize 
the main model assumptions, with a focus on several 
modifications. 

2.1. Equations and Assumptions 

In our model, CRs and hot gas are treated as an in- 
teracting "two-fluid" system, whose dynamical evolution 
may be described by the following four equations: 



dp 
'dt 



pV • V = 0, 



(1) 



= - V(P + Pc) - pV^ + V • n, (2) 



de 

— + V • (ev) = -PV • V + n : Vv, (3) 



dec 
dt 



+ V • (ecv) = -PeV • V + V • (/^Vec), (4) 



where d/dt = d/dt-\-v -S/ is the Lagrangian time deriva- 
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tive, K, is the CR diffusion coefficient, H is the viscous 
stress tensor (see the Appendix), and ah other variables 
have their usual meanings. The gas pressure P and CR 
pressure Pc are related with the gas and CR energy den- 
sity e and Cc via P = (7 — l)e and Pc = (7c — l)ec 
respectively, where we assume 7 = 5/3 and 7c = 4/3. 
The nature of the relativistic particles with energy den- 
sity Cc is unspecified and may be electrons and/or pro- 
tons with any spectra. Of course the equation of state 
may be somewhat harder if Cc is mainly contributed by 
trans-relativistic protons at ~ 1 GeV. 

We assume a temporally fixed Galactic potential <l>, 
which is contributed by three components: the bulge, 
disk and dark matter halo. Their properties are elabo- 
rated in Section 2.3 of Paper I. At time t = we as- 
sume that the CR energy density is zero in the Galaxy, 
ec = 0. The hot Galactic gas is assumed to be initially 
isothermal with temperature T = 2.4 x 10^ K and its 
density distribution is solved from the assumption of hy- 
drostatic equilibrium. The free parameter ngo, the ther- 
mal electron number density at the origin, determines 
the normalization of the gas density distribution and is 
chosen to be 0.1 cm~^ throughout this paper. The re- 
sulting initial density distribution is shown in the top 
panel of Figure [TJ while the current density distribution 
after the AGN event in a typical run V3-diff3 is shown 
in the bottom panel. The dependence of our results on 
the poorly-constrained parameter Uqo has been explicitly 
explored in Paper I, which shows that for higher values of 
neo, more powerful jets are needed to produce the Fermi 
bubbles with the same observed morphology. 

We ignore radiative cooling of thermal gas, which is 
unimportant during the short-duration (< 1-3 Myr) of 
our simulations. We also neglect CR energy losses from 
synchrotron and IC emissions, which is a good assump- 
tion for CR protons and CR electrons at 10 - 100 GeV or 
lower energies. In this paper, we follow the evolution of 
the integrated CR energy density Cc (equation |4|) , which 
may not be significantly affected by the CR cooling as it 
is probably mainly contributed by low-energy CR elec- 
trons and possibly CR protons with long lifetimes. The 
cooling may be important for electrons at TeV energies, 
and thus may significantly affect the gamma-ray spec- 
trum if TeV electrons dominate gamma-ray emissions 
from the Fermi bubbles (in particular at high latitudes). 
The main purpose of this paper is to reproduce the basic 
observed morphology of the gamma-ray bubbles, but we 
defer a detailed direct comparison with the data (includ- 
ing spectral predictions) to future work. 

Equation |4] describes the evolution of CR energy den- 
sity including both advection and diffusion. Note that 
- as discussed in Section 2.1 of Paper I - we ignore 
CR streaming, which may play a similar role as CR 
diffusion, transporting CRs away from local thermal 
plasma. Typical values of CR diffusivity hi in the 
Galaxy are found to be hc ^ (3 — 5) x 10^^ cm^ s~^ 
for CRs at about 1 GeV (jStrong et al.l I2007D . How- 
ever, in Paper I we show that to produce the sharp 
edges of the observed Fermi bubbles, the CR diffusivity 
across the bubble surface must be strongly suppressed, 
which may occur naturally if magnetic field lines are 
mainly tangential on the bubble surface, as expected 
from magnetic draping. Previous work has explored 
similar effects in AGN blown bubbles in galaxy clus- 



ters, with similar results (jMathews fc Brighentil 120071 : 
iRuszkowski et all 12008). Thus for most runs, we choose 
a strongly-suppressed uniform and constant CR diffusiv- 
ity = 3 x 10^^ cm^ s~^, which ensures that CR diffu- 
sion has a negligible effect on the few-Myr evolution of 
our simulated bubbles. In Section 3.2, we present three 
additional runs where without increasing CR diffusivity 
outside the bubble, we increase CR diffusivity in the bub- 
ble interior to tvi^t = (1 — 6) x 10^^ cm^ s~^, exploring 
the effect of CR diffusion on the CR distribution in the 
bubble interior. 

Equations ([I]) — (|4]) are solved in (i?, z) cylindrical co- 
ordinates usi ng a 2D axisymmetric E ulerian code similar 
to ZEUS 2D (jStone Normanlll992[ ).'^ In particular, we 
have implemented important new physics into the code, 
including gas viscosity (see the Appendix), CR advection 
and diffusion, and the dynamical interaction between hot 
gas and CRs. The computational grid consists of 400 
equally spaced zones in both coordinates out to 20 kpc 
plus additional 100 logarithmically-spaced zones out to 
50 kpc. The jet inflow is introduced along the 2;-axis (the 
rotation axis of the Galaxy) from the GC. See Paper I for 
the details of our jet injection method. The jet parame- 
ters are the same in all the runs presented in this paper: 
speed Vjet = 3.0 x 10^ cm/s, radius Rjet = 0-4 kpc, dura- 
tion tjet = 0.4 Myr, CR energy density ejcr = 1.0 x 10~^^ 
erg/cm^, thermal gas density pj = 1.102 x 10~^^ g/cm^ 
(density contrast r] = 0.01 with respect to the ambient 
gas), thermal energy density ej = 5.4 x 10~^^ erg/cm^. 
These jet parameters produce CR bubbles having similar 
morphologies as the observed Fermi bubbles. The total 
power of this jet is Pjet ~ 8.6 x 10^^ erg s~^, dominated 
by the kinetic power. The total energy injected by these 
two opposing jets is 2£^jet = 2Pjet^jet ^ 2.17 x 10^^ erg. 
We stop each simulation at time t = ^Fermi, when the 
produced CR bubble reaches z = 10.5 kpc along the jet 
axis. At this time, the projected CR bubble viewed in 
the Galactic coordinate system roughly reaches the high- 
est latitude to which the observed Fermi bubbles extend 
(see the bottom panels of Fig. [6]). Thus ^Fermi is the 
predicted dynamical age of the Fermi bubbles in each 
model, and depends on model parameters. 

The dependence of our results on jet parameters has 
been explored in Paper I. In particular, we point out that 
the CR energy density ejcr is roughly degenerate with 
thermal energy density ej in the sense that the jet evolu- 
tion depends on the total jet pressure (ejcr/3 + 2ej/3) and 
is not very sensitive to the specific CR pressure. Thus 
our current model can not uniquely predict the gamma- 
ray luminosity of the resulting Fermi bubbles. However, 
the observed gamma ray flux may be used to constrain 
the particle content in the bubbles. We did some sim- 
ple emission calculations in Section 3.5 of Paper I and 
found that (1) If the gamma-ray emission from the bub- 
bles is mainly produced by CR electrons, the required CR 
electron pressure is negligible compared to the total bub- 
ble pressure, which may instead be dominated by other 

^ Since both vorticity and magnetic flux are subject to stretching 
in 3D, but not in 2D, the restriction to 2D has more of an effect in 
MHD simulations of the KH instabihty, which attempt to charac- 
terize th e fleld strength r equired for magnetic tension to stabilize 
the flow (|Rvu et al.ll200Cll ). an issue beyond the scope of this study. 
Such studies generally flnd that a smaller initial fleld is required in 
3D to stabilize the flow, due to greater fleld ampliflcation. 
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TABLE 1 
List of Simulations 



Mvisc ^int ^Fermi 

Run (g cm~^ s~^) (cm^ s~^) (Myr) 

VO 3 X 10^^ 1.85 

V0.5 0.5 3 X 10^6 1.81 

VI 1 3 X 10^6 1.87 

V3 3 3 X 10^6 1.67 

VIO 10 3 X 10^6 1.60 

V30 30 3 X 1026 2.24 

V3-diffl 3 1 X 1028 1.60 

V3-diff3 3 3 X lO^s 1.49 

V3-difr6 3 6 X lO^s 1.39 



^^int is the value of CR diffusion coefficient within the evolv- 
ing CR bubble (see Sec. 3.2 for details). The CR diffusivity 
outside the bubble is always chosen to be 3 x 10^^ cm^ s"-*^, an 
arbitrarily-chosen small value to suppress CR diffusion across 
bubble edges. 

components, e.g., thermal gas, CR protons, or magnetic 
fields. (2) If the gamma-ray emission is mainly due to 
CR protons, the required CR proton pressure is much 
higher, probably dominating the total bubble pressure. 

2.2. The Role of Viscosity in Jet/Bubble Evolution 

One of the main goals of this paper is to study the role 
of shear viscosity on the evolution of the Fermi bubbles 
in the jet scenario. This study is mainly motivated by our 
previous jet simulations presented in Paper I, where the 
Kelvin-Helmholtz (KH) and potential Ray leigh- Taylor 
(RT) instabilities develop at the edges of the resulting 
CR bubbles, strikingly inconsistent with smooth edges 
of the observed Fermi bubbles. Smooth bubble edges 
suggest that the instabilities are effectively suppressed 
by some additional physics, which works on small scales 
along the whole bubble surface. Viscosity is an ideal 
candidate mechanism to fulfill this purpose. 

In a fully ionized, unma^neti zed, thermal plasma, the 
dynamic viscosity coefficient is (|Braginskiilll958l : iSpitzerl 
[1962.): 

/I A\~^ / T \ 
Mvisc=6.0xl03(^^j ^j^j gcm-is-l:?) 

where T is the temperature of the plasma in Kelvin and 
In A is the Coulomb logarithm. One important property 
of the viscosity is that it increases dramatically with gas 
temperature (/ivisc ^ T^/^). For example, the value of 
/ivisc increases from 0.06 to 6000 g cm~^ s~^ when tem- 
perature increases from 10^ to 10^ K. As shown in Paper 
I, the AGN jet event induces a strong shock propagating 
into the hot halo gas, which heats the gas to tens to hun- 
dreds of keV at early times. The gas temperature drops 
as the gas expands into the halo, but even at the current 
time, the shocked gas still has temperatures of a few keV. 
Large viscosity in such hot gas may potentially play a sig- 
nificant role in the evolution of the Fermi bubbles (more 
generally AGN bubbles), in particular, suppressing the 
development of KH and RT instabilities. Unknown mag- 
netic fields may suppress viscosity across field lines, but 
as we show in this paper, a very small fraction (less than 
1%) of the viscosity in equation (|5|) is capable of sup- 
pressing these instabilities. 
Unlike the theories of accretion disks, where the role 



of viscosity is greatly appreciated, theoretical/numerical 
studies of AGN jets often ignore viscosity. One rea- 
son for this neglect is the difficulty in attributing ob- 
servational features directly with the effects of viscos- 
ity; the smooth edges of the Fermi bubbles may provide 
an unusual opportunity to place such observational con- 
straints. Similar smooth edges have also been observed 
in m any radio bubbles in g alaxy clusters, which moti- 
vated iRevnolds et al.l (|2005h to study the role of viscos- 
ity on the evolution of buoyantly rising bubbles. These 
studies also found that modest levels of viscosity could 
stabilize Ray leigh- Taylor and Kelvin-Helmholz instabil- 
ities, allowing the bubbles to m aintain their in t egrity . 
However, the numerical study of IRevnolds et al.l (|2005l ) 
only considered the buoyant rise of initially static bub- 
bles, side-stepping the initial jet-driven inflation of the 
bubble. As they acknowledge, besides excluding the ef- 
fect of the jet on its surroundings (such as driving strong 
shocks), this leaves out complex internal motions within 
the bubble which arise during the inflation phase. By 
directly simulating the AGN jet, we take such effects 
into account. We find that internal backflows within the 
bubble contribute strongly to the development of fluid in- 
stabilities at the bubble surface, and cannot be ignored; 
indeed, viscosity plays a critical role in mitigating such 
backflows. 

In this paper, we adopt a constant, isotropic viscos- 
ity and study how the results vary with different val- 
ues of viscosity in a series of simulations (see Table 1). 
This simplifled appr oach, which was also adopted by 
IRevnolds et al.l [20051 allows us to make a preliminary 
assessment of how bubble evolution might be affected 
by viscosity (see Figure [2j). The true nature of viscos- 
ity here is highly uncertain. For instance, equation (j5j) 
describes the isotropic viscosity coefficient in an unmag- 
netized hot plasma, but the gas/plasma in and outside 
the Fermi bubbles contains magnetic fields, which makes 
viscosity anisotropic, as it is enormously suppressed (by 
a factor ~ 10^^) across field lines. The exact value of vis- 
cosity in a turbulent medium with tangled field lines is 
unknow n, although perhaps in analog y with thermal con- 
duction ([Naravan fc MedvedevI 120011 ). values - 1 - 30% 
of the Braginskii-Spitzer value are plausible. The value 
could be considerably smaller if the field is coherent: for 
instance, nearly-parallel magnetic field lines at the bub- 
ble surface, as suggested by sharp bubble edges (see Sec. 
3.2 of Paper I), may significantly suppress momentum 
transport across bubble edges, although some tangling 
of the field here-perhaps due to the instabilities itself- 
could still allow a non- negligible value. The nature and 
level of viscosity in the bubble interior are even more un- 
certain, as the thermal gas there is extremely hot and 
underdense; the formal Coulomb mean free path: 

is so large that it is effectively collisionless. 

Thus, while we compare the values of viscosity that 
we use to the Braginskii-Spitzer value to illustrate its 
magnitude, we make no claims as to its provenance. For 
instance, the stress could be magnetic in nature. It could 
also arise from anisotropic pressure, which is sensitive not 
just to the topology but the magnitude of magnetic fields. 
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Fig. 2. — Central slices (16 x 15 kpc) of CR energy density in logarithmic scale in runs VO, V0d5, VI, V3, VIO, and V30 at t = ^Fermi? 
which is shown at the top of each panel for the corresponding run. Horizontal and vertical axes refer to R and z respectively, labeled 
in kpc. The stabilizing effect of viscosity on bubble edges can be clearly seen here as viscosity increases from panel to panel, and the 
Kelvin-Helmholtz and Ray leigh- Taylor instabilities are fully suppressed when //vise ^ 3 g cm~^ s~^. 



In a weakly collisional/collisionless plasma such as the 
bubble interior, pressure anisotropy p\\ ^ arises from 
conservation of the magnetic moment for each particle 
II — mv\/2B = const, which implies that any change in 
the field is accompanied by a change in the perpendicular 
pressure to keep p^/B '-«^const. This then triggers micro- 
instabilities (such as the firehose, mirror, ion cyclotron 
instabilities) which feed off the press ure anisotropy an d 
pin it at marginal stability values (|Rosin et alJ l2Qlll ). 
The micro-instabilities change the pressure anisotropy 
either via an enhanced rate of collisions through an ef- 
fective particle scattering mechanis m, a source of effec- 
tive viscosity (|Sharma et alJ [2'QQ6f ), or modification of 
the rate of strain of the magnetic field so as to cancel 
the pressure anisotropy created by the cha nging fields 
(jRosin et alJl2Qlll : ISchekochihin et all 1201 Q[ ): the latter 
gives rise to a viscosity in a turbulent medium that scales 
as the parallel Braginskii value (and by dissipating tur- 
bulent motions, c ould provide significant viscous heating; 
iKunz et alJl2Qlll ). Viscosity in collisionless plasma may 
also be caused by particle scattering with magnetic ir- 
regularities and Alfven waves, which has been invoked 
to explain the origin of CR diffusion - a well-known 
transport process in collisionless plasma. Assuming that 
Mvisc ^ pvX^ the effective mean free path of proton scat- 



tering for our assumed level of viscosity is: 



A ~ 1 kpc 



3 g cm ^ s 



'—){ 



. 10^ cm s" 



lQ-29 g cm~3 



(7) 



where v is the kinetic velocity of protons and p is the 
plasma density. 

Thus, while the nature of viscosity in this context is 
highly uncertain, assuming an isotropic, uniform vis- 
cosity is not unreasonable. The next step would obvi- 
ously be to perform MHD simulations similar to those of 
(|Sh arma et aIll2QQ6D for accretion disks. It would be ex- 
citing to place empirical constraints on viscosity based on 
comparisons of our calculations with the observed Fermi 
bubbles. 

In the Appendix, we explicitly present our numerical 
method to implement the fully compressible shear vis- 
cosity into our 2D code. The viscous runs are fairly 
expensive, because the time-step imposed by viscosity 
scales with p(Ax)^//ivisc7 where Ax is the resolution of 
the computational grid. In particular, the viscous time- 
step becomes extremely small at some small regions in 
the bubble interior, where the thermal gas density is very 
low due to the low initial jet density, the bubble ex- 
pansion and viscous heating. To allow the simulations 
to proceed, we thus turn off viscosity in computational 
cells where the thermal gas density drops below 10~^^ 
g cm~^. This restriction only affects some small regions 
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deep inside the bubbles, and does not appreciably affect 
the bubble evolution. 

3. RESULTS 

Here we present the results of our numerical calcula- 
tions, which are compared directly with the gamma-ray 
observations of the Fermi bubbles. Our main purpose is 
not to perfectly reproduce the Fermi bubbles, but rather 
we aim to identity potential physical mechanisms rele- 
vant for the Fermi bubble event. In particular, we inves- 
tigate the potential roles of viscosity and CR diffusion 
on the bubble evolution and the CR distribution within 
the bubbles. 

3.1. Suppression of KH Instabilities 

To study if viscosity can indeed suppress the develop- 
ment of KH instabilities, we performed a series of simu- 
lations with different levels of viscosity. These runs are 
denoted as "run V/i", where the viscosity coefficient for 
our six runs are chosen to be /ivisc — 

0, 0.5, 1, 3, 10, and 
30 g cm~^ s~^ (see Table 1). In these runs, CR diffusion 
plays a negligible role, as CR diffusivity is always chosen 
to be a small constant = 3 x 10^^ cm^ s~^ (see equation 
7). 

Figure [2] shows central slices of CR energy density dis- 
tribution in these runs at t = ^Fermi , when the resulting 
CR bubble reaches a distance of 2: = 10.5 kpc along the 
jet direction. In the non- viscous run VO, irregularities 
clearly develop at the surface of the jet-induced CR bub- 
ble, indicating the significance of KH instabilities. As the 
viscosity coefficient increases in runs V0.5 and VI, the 
magnitude of surface irregularities becomes smaller. In 
run V3, VI 0, and V30, surface irregularities disappear. 
Thus viscosity effectively suppresses RT and KH insta- 
bilities when //vise ^ 3 g cm~^ s~^. This is exactly what 
we expected: As the viscosity coefficient increases, the 
viscous stress at the bubble surface becomes more signif- 
icant, effectively suppressing the growth of small pertur- 
bations, which would otherwise grow into large vortices. 

As discussed in Section 2.2, the true value of viscosity 
in and out the Fermi bubbles is very unclear, particularly 
due to the uncertain role of magnetic fields in transport 
processes. The values of viscosity adopted in success- 
ful runs V3, VIO, and V30, respectively /ivisc — 

3, 10, 

30 g cm~^ s~^, are much less than the Spitzer viscos- 
ity of hot gas surrounding the bubbles. This hot gas 
is strongly heated by the jet-induced shock, and cur- 
rently has temperatures of ~ 0.6 - 4 x 10^ K (even 
higher at earlier times). For comparison, the temper- 
ature 10^ K corresponds to the Spitzer viscosity (eq. 
5) of ~ 6000 g cm~^ s~^. Thus the adopted viscosities 
in these runs are only around 0.1% - 1% of the Spitzer 
value in the surrounding gas. Such a low viscosity level 
may represent the true level of momentum transport rate 
across bubble edges. This is due to the suppression of 
transport processes across nearly-parallel magnetic field 
lines at the bubble surface, as discussed in Section 2.2. It 
is likely that near the bubble surface, the magnetic field 
lines are not completely tangential and a small level of 
field tangling allows a low level of momentum transport 
across the bubble surface, which is strong enough to sup- 
press the instabilities as shown in our calculations. The 
importance of viscous transport inside but near the bub- 
ble surface can be seen by small values of the Reynolds 



number in the jet backflow: 



Re -0.2 



10-29 g cm-^J V2OOO km s"^ 

1 



L 



1 ,,-1 



0.1 kpc y \3 g cm ^ s 



where p ^ 10 g cm ^ and L 



(8) 



0.1 kpc are roughly 
the thermal gas density and thickness of the jet backflow 
layer. In contrast, the Reynolds number in the ambient 
shock- heated gas is much larger: 



Pvis 



L 



-3 

Pvisc 



10-28 g J \2000 km/s 
1 



10 kpc J \3 g cm~^ s~^ 



(9) 



where p — 10~^^ g cm~^ and L ~ 10 kpc are roughly 
the thermal gas density and velocity length scale in the 
shocked halo gas, indicating that viscosity at this level is 
dynamically unimportant there. 

Viscosity also helps dissipate gas motions in the bub- 
ble interior, reducing the level of CR advection, as clearly 
seen in Figure [21 When viscosity is not important (e.g., 
in runs VO and V0.5), CR advection driven by circulating 
gas motions transports and mixes CRs within the bub- 
ble, producing a volume-filling CR bubble. In contrast, 
when viscosity becomes important (in runs V3, VIO and 
V30), gas motions are significantly reduced, and CRs are 
mainly located in the jet backflow near bubble edges. In 
the runs with low viscosity, the jet backflow forms a large 
scale circulating flow in the bubble interior, as better seen 
in run VI (the right-top panel), which advect CRs, pro- 
ducing a volume-filling CR bubble. However, in runs V3, 
VIO, and 30, viscosity dissipates the shear motions in- 
duced by the backflow near bubble edges, preventing the 
formation of circulating motions in the bubble interior. 

This viscosity effect can be better seen in Figure [H 
which compares the early-stage jet evolution in runs VO 
and V3. Here the velocity field is over-plotted on the 
distribution of CR energy density. As clearly seen, in 
the non- viscous run VO, the jet backflow induces signifi- 
cant circulating motions, which advect CRs to the bubble 
interior. In contrast, in the viscous run V3, circulating 
motions are only seen at the very early time t = 0.1 Myr, 
and at later times, the jet backflow stops flowing back- 
ward and mainly expands into the Galactic halo. One 
consequence of this viscosity effect is the reduction of 
shear motions at the bubble surface, which helps sup- 
press the growth of KH instabilities. 

3.2. The Flat Gamma-ray Surface Brightness 

An important feature of the observed Fermi bubbles 
is the approximately uniform gamma-ray surfa ce bright- 
ness, particularly at high latitude (|&| ^ 30°; see lSu et al.l 
I2OIOI ). It is not trivial to form such a flat surface bright- 
ness distribution. The top panel of Figure |4] shows the 
line-of-sight projected CR energy density distribution 
in Galactic coordinates in the non-viscous run VO at 
t = tpermi- I^ addition to instability-induced irregular- 
ities at bubble edges, the bubble is clearly seen to be 
center-brightened. Indeed, a spatially-uniform CR dis- 
tribution produces a center-brightened surface brightness 
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Run VO, t = 0.1 Myr 
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Run V3, t = 0.1 Myr 




Run VO, t = 0.3 Myr 



Run V3, t = 0.3 Myr 







Run VO, t = 0.5 Myr 




Fig. 3. — Central slices of CR energy density in logarithmic scale in run VO (top panels), and V3 (bottom panels) at t = 0.1 (left), 0.3 
(middle) and 0.5 Myr (right). Horizontal and vertical axes refer to R and z respectively, labeled in kpc. Arrows superposed show gas 
velocity. White dots within the CR bubble indicate that the local velocity exceeds 2.5 x 10^ cm/s, while those far away from the bubble 
indicate vanishing velocities. In the non-viscous run VO, the jet backflow at the bubble edges reaches the bottom, and then rises up in the 
interior, forming a large scale circulation. However, in the viscous run V3, the jet backflow stops flowing backward quickly and stays at 
the expanding surface without forming circulations. 



in projection. The roughly flat surface brightness implies 
that the CR density gradually increases toward the bub- 
ble edge, and particularly increases with latitude from 
the bubble center. 

It is of great interest to study why the CRs should 
concentrate near the bubble edges. CR diffusion and 
streaming tend to make the CR distribution more uni- 
form, and CR advection associated with circulating mo- 
tions in the non- vi scous run VO can also ni ix CRs deep 
within the bubble. iMertsch fc Sark eir (2011) argued that 
stochastic re-acceleration of CR electrons could account 
for the flat gamma-ray surface brightness if CR electrons 
are preferentially re- accelerated near bubble edges. They 
argued that KH instabilities seen in our non- viscous run 
VO produce turbulence, resulting in CR re-acceleration 
preferentially near bubble edges. However, it is impor- 



tant to note that the instabilities seen in our non- viscous 
runs represent large scale turbulence while the 2nd-order 
Fermi acceleration mechanism is dominated by small 
scale turbulence. Furthermore, generically this mecha- 
nism also produces limb brightened profiles, and the dif- 
fusion properties of the bubble interior must be adjusted 
to fit the data. Ultimately though, the Fermi observa- 
tions show that the bubble edges are very smooth, sug- 
gesting_that the level of turbulence is low in these regions 
(jSu et al.l[20Toh . 

3.2.1. The Potential Role of Viscosity 

In the previous subsection, we show that a small level 
of viscosity (compared to the Spitzer viscosity in the sur- 
rounding gas) can suppress the development of KH insta- 
bilities, resulting in smooth bubble edges. Furthermore, 
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og(projected e^), run VO, t=1.85 Myr 




bg(projected e^,), run V3, t=1.67 Myr 

80 60 40 20 -20 -40 -60 -80 

8.50 8.92 9.33 9.75 10.17 10.58 11.00 

Fig. 4. — Line-of-sight projected CR energy density in logarith- 
mic scale in run VO (top) and V3 (bottom) at t = tpermi- Hor- 
izontal and vertical axes refer to Galactic longitude and latitude 
respectively, labeled in degrees. The dotted region in each panel 
encloses the observed north Fermi bubble, while the solid circle en- 
closes the south Fermi bubble. Edge irregularities are clearly seen 
in the non-viscous run VO, while the observed Fermi bubbles show 
smooth edges. The viscous run V3 shows smooth edges, but the 
gamma-ray intensity distribution is limb-brightened, inconsistent 
with the observed flat surface brightness. 

viscosity reduces the levels of gas shear motions and the 
associated CR advection in the bubble interior, signifi- 
cantly affecting the spatial distributions of CRs and the 
line-of-sight projected gamma-ray intensity in the Fermi 
bubbles. Due to its relatively small inertia, the jet is 
deflected at the top of the bubble and flows backward, 
transporting CRs down along the bubble boundary. In 
the absence of viscosity this backflow is deflected once 
again at the bubble bottom and returns upward in the 
direction of the original jet, filling most of the bubble 
interior with CRs. However, in viscous runs (e.g. run 
V3) this second upward motion is damped by viscosity 
and the original boundary backflow is arrested and ex- 
pands with the bubble, retaining the concentration of 
CRs near the bubble boundary produced by the original 
backflow as seen in the bottom panels of Fig. [2j Thus, 
an edge-favored CR distribution, inferred from the ob- 
served flat gamma-ray intensity, is a natural consequence 
of shear viscosity, which is often ignored in previous jet 
studies. An edge-favored CR distribution may also be 
present in some extragalactic radio bubbles, where the 
observed r adio synchrotron emissivity is peaked a t bub- 
ble edges (jCarvalho et al.l I2QQ5I : iDalv et al.llMTol) . We 
note that it is difficult to explain this observational fea- 
ture by other physical mechanisms. 

Momentum transport near the bubble surface is critical 
in suppressing the backward motion of the jet backflow. 
Figure [5] shows variations of Vz and CR pressure along 
the i?-direction for the non-viscous run VO (left panels) 
and the viscous run V3 (right panels) at z = 2 kpc at 
three times t = 0.3, 0.5, and 1.0 Myr. The jet backflow 
layer is located at the bubble surface, corresponding to 



the CR pressure peak at the right end of each line in bot- 
tom panels of Figure [51 In top panels, it is represented 
by regions with negative values of Vz . Due to the bubble 
expansion, the backflow moves to larger Galactocentric 
distances with time. Comparing the backflow layers in 
the left and right panels, it is clear that viscosity signif- 
icantly suppresses the backward motion of the backflow 
in run V3. 

The top panels of Figure [5] clearly show the presence of 
velocity gradients at both the inner and outer surfaces of 
the backflow, suggesting that momentum transport into 
the backflow from both the bubble interior and ambient 
gas contributes to the suppression of its backward mo- 
tion in run V3, which adopts a spatially constant viscos- 
ity coefficient. It is possible that momentum transport 
across one surface alone may be sufficient to reduce the 
backflow's backward motion. In particular, strong veloc- 
ity gradients are present near the inner interface of the 
backflow with the bubble interior as clearly seen in both 
the left-top (non-viscous) and right-top (viscous) panels. 
We speculate that momentum transport across the in- 
ner interface alone is sufficient to reduce the backflow's 
backward motion and thus suppress KH instabilities if 
momentum transport across the outer interface is fully 
suppressed by parallel magnetic fields. We tentatively 
confirm this speculation in a few additional simulations 
where viscosity is only allowed in the bubble interior. 
However, it is not easy in these simulations to accurately 
determine the bubble surface (i.e., to fully shut off mo- 
mentum transport across the bubble surface), which is 
resolved by a few numerical cells. More robust conclu- 
sions can only be made by future high resolution simu- 
lations with more advanced numerical technologies. 

3.2.2. The Potential Role of CR Diffusion 

The flat gamma-ray surface brightness is a very in- 
triguing observational feature, requiring a gradual in- 
crease of CRs toward the bubble surface. As noted pre- 
viously, a uniform CR distribution will give rise to a 
center-brightened gamma-ray surface brightness. If most 
CRs are concentrated at the bubble surface, the gamma- 
ray surface brightness will instead be limb-brightened, 
as clearly seen in the bottom panel of Figure H] which 
shows the distribution of line-of-sight projected CR en- 
ergy density in the Galactic coordinate system in a typ- 
ical viscous run V3 at t = tFermi- We now show that 
including CR diffusion within the bubbles can lead to 
a profile which roughly matches observations. In run 
V3, we choose a spatially uniform low CR diffusivity 
{k, = 3 X 10^^ cm^ s~^), but as shown in Paper I, CR 
diffusion only needs be suppressed across the bubble 
surface to reproduce the observed sharpness of bubble 
edges. Such suppression could be produced by mag- 
netic draping (particularly at early times), which pro- 
duces magnetic field lines appr o ximately tangent to the 
bubble surface fLvutikov| 120061 : iRuszkowski et al.l 120071 : 
iDursi fc Pfromm er 200l). No such considerations ap- 
ply to the bubble interior, whose field structure is likely 
set by the bubble inf lation process and subsequent re- 
connection (e.g., see iBraithwaitd [2010f ). probably not 
significantly suppressing CR diffusion there. We note 
that the details of the potential magnetic draping asso- 
ciated with the supersonic GC jets here may be differ- 
ent from those described in IRuszkowski et al.. (^2007i) and 
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Fig. 5. — Variations of the z-component gas velocity (top) and CR pressure (bottom) along the i?-direction (perpendicular to the jet 
axis) for the non-viscous run VO (left panels) and the viscous run V3 (right panels) dX z = 2 kpc at three times t = 0.3, 0.5, and 1.0 Myr. 
The jet backflow is represented by regions with negative values oi Vz in top panels, corresponding to CR pressure peaks at the right end 
of each line in the bottom panel. The backflow layer is located at the bubble surface and follows the global expansion of the CR bubble. 
In run V3, its backward motion is reduced by momentum kinetically transported from both the bubble interior and the ambient halo gas. 
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Fig. 6. Top: central slices (16 x 15 kpc) of log(ec) in run V3-diffl (left), V3-diff3 (middle), and V3-diff6 (right) at t = trermi- Horizontal 
and vertical axes refer to R and z respectively, labeled in kpc. Arrows superposed show thermal gas velocity. Bottom: line-of-sight projected 
CR energy density in these three runs at t = tpermi- Horizontal and vertical axes refer to Galactic longitude and latitude respectively, 
labeled in degree. The CR diffusivity in the bubble interior, which increases from left to right, may play an important role in explaining 
the flat gamma-ray surface brightness of the Fermi bubbles. 
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iDursi fc Pfrommerl (|2QQ8f ), where the draping is caused 
by subsonic bubble motions. Clearly, MHD jet simula- 
tions, which are beyond the scope of the current paper, 
would be very helpful in understanding the evolution of 
the magnetic field topology near the bubble surface. 

To accurately study CR diffusion during the jet/bubble 
evolution, we need to rely upon future magnetohydrody- 
namic (MHD) simulations with anisotropic CR diffusion. 
However, within our current methodology, we can still 
explore the potential role of CR diffusion on the CR dis- 
tribution in the bubble interior. To this end, here we 
present three additional runs V3-diffl, V3-diff3, and V3- 
diff6, in which we increase the value of CR diffusivity in 
the bubble interior to /^int = 1 x 10^^, 3 x 10^^ and 6 x 10^^ 
cm^ s~^ respectively (see Table 1 for other model pa- 
rameters). We distinguish the bubble interior (p < Pcrit) 
from the outside regions (p > Pcrit) surrounding the ex- 
panding bubble by a density criterion, as the CR bubble 
is separated from the surrounding halo gas through a 
contact discontinuity, across which thermal gas density 
increases abruptly (see the bottom panel of Fig. [T]). We 
use pcrit to identify gas inside the bubbles (low p, high k,) 
from ambient shocked gas (high p, low tv). Since in our 
model there are essentially no CRs outside the bubbles, 
the low CR diffusivity there only suppresses CR diffu- 
sion across the bubble surface. As the bubble expands 
quickly, the thermal gas density within the bubble drops. 
Therefore, the critical density (pcrit) identifying the bub- 
ble interior should also drop with time. We choose pcrit to 
be twice the volume- averaged thermal gas density along 
the jet axis within the CR bubble at any time during the 
bubble evolution. Such a crude method produces accept- 
able results, as clearly seen in the top panels of Figure 
[6l which shows spatial distributions of CR energy den- 
sity at t = tpermi in ruus V3-diffl, V3-diff3, and V3-diff6. 
The edges of the resulting bubbles are as sharp as in the 
low- diffusivity run V3, and the bubble morphology is also 
similar in all these runs. This is consistent with what we 
would expect, since CR diffusion is only increased in the 
bubble interior while still significantly suppressed across 
the bubble surface. As the density jump across the bub- 
ble surface is quite large, our results are not sensitive to 
the specific value of pcrit- 

CR diffusion transports CRs near bubble edges to the 
bubble interior, particularly during early times when the 
size of the bubble is small. The typical length that CRs 
diffuse within a duration of t is / ^ x/kI: 



I - 0.3 



3 X 1028 cmVs 



1/2 



1 Myr 



1/2 



kpc. (10) 



As A^int increases, the CR distribution within the Fermi 
bubbles becomes less limb-brightened and more uniform, 
as clearly seen in the top panels of Figure [6] (from left to 
right). The bottom panels of Figure [6] show the line-of- 
sight projected CR energy density in runs V3-diffl, V3- 
diff3, and V3-diff6 in Galactic coordinates at t = tpermi- 
In run V3-diffl with tZi^it = 1 x 10^^ cm^ s~^, the dis- 
tribution of projected CR energy density is still limb- 
brightened, similar to the low- diffusivity run V3. But 
in runs V3-diff3 and V3-diff6 with /^int = 3 x 10^^ and 
6 X 10^^ cm^ s~^ respectively, the projected CR energy 
density distribution becomes very flat at high latitude 
(|6| > 30°), consistent with the gamma-ray observations 




-10 -20 

longitude 

Fig. 7. — Line-of-sight projected CR energy density in run V3- 
diff3 at t = tpermi cis a function of Galactic longitude (in degrees) 
at three latitudes: 6 = 20°, 30°, 40°. The projected CR distribu- 
tion is roughly uniform with longitude, and slightly increases with 
latitude. 

of the Fermi bubbles. This can also be seen in Figure 
[7J which shows longitudinal variations of the line-of-sight 
projected CR energy density in run V3-diff3 at t = tpermi 
at three latitudes: 6 = 20°, 30°, 40°. At lower latitude, 
the projected CR energy density is slightly lower, but 
higher ISRF intensities there could boost the gamma-ray 
IC emissivity. The flatness of the gamma ray intensity 
with latitude, which needs to be further corroborated by 
three-year or longer Fermi observations, seems to require 
a fine-tuning of the latitudinal distribution of CR parti- 
cles. The effects of non-uniform ISRF on the projected 
gamma-ray intensity will be explored in future work. 

In the discussions above, we have mainly considered 
gamma ray emissions due to CR electrons through IC 
scattering. If CR protons are also present in the Fermi 
bubbles, they will also produce gamma ray emissions. 
It is yet unclear if the gamma ray emission from the 
Fermi bubbles is do minated by CR electrons or p rotons 
(|Dobler et all 120101 : iCrocker AharonianI [201 If ). The 
line-of-sight projections of pCc at t = tpermi in runs V3 
and V3-diff3 shown in Figure [8] indicate that in these vis- 
cous runs (also seen in runs V3-diffl and V3-diff6) the 
gamma-ray surface brightness contributed by CR pro- 
tons peaks at the bubble edge, where the jet backflow 
is located. This edge concentration occurs because the 
thermal gas density in the jet backflow is much higher 
than that in the expanding bubble interior, as seen in 
the bottom panel of Figure [H The gas density in the 
jet backflow is mainly determined by the initial jet den- 
sity, which is not well constrained in our current model 
(see Section 3.3 and 3.4 of Paper I). But some level of 
fine-tuning of the initial jet density may be required if the 
relatively flat gamma-ray surface brightness is dominated 
by CR protons in our viscous jet scenario. Further stud- 
ies are required to investigate if the gamma ray emission 
of the Fermi bubbles is dominated by CR electrons or 
protons, which is beyond the scope of the current paper. 

3.3. The Evolution of the Fermi Bubbles 

In this subsection, we briefly comment on the temporal 
evolution of the Fermi bubbles in our model. While CR 
diffusion transports CRs to the bubble interior, leading 
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Fig. 8. — Line-of-sight projection of pec (a proxy for the gamma- 
ray surface brightness originated from CR protons) in logarithmic 
scale in run V3 (top panel) and V3-diff3 (bottom) at t = ^Fermi- 
Horizontal and vertical axes refer to Galactic longitude and latitude 
respectively, labeled in degrees. The gamma-ray surface brightness 
contributed by CR protons peaks at the jet backflow near the bub- 
ble surface, where thermal gas density is much higher than that in 
the bubble interior, as seen in the bottom panel of Figure [T] 

to a spatially-uniform CR distribution, CR advection as- 
sociated with the fast bubble expansion transports CRs 
outward, tending to counteract this effect. Thus, a CR 
distribution which increases toward the bubble surface, 
producing a roughly flat gamma-ray surface brightness, 
exists for a quite long time (at least ~ 1 Myr) during 
the bubble evolution. Figure [9] shows CR distributions 
(top) and the line-of-sight projected CR distributions in 
run V3-diff3 at times t = 0.7, 1, and 2 Myr. Clearly at 
t = 0.7 Myr, the CR concentration at the compressed jet 
tip can still be seen, but at time t = 1 - 2 Myr, the pro- 
jected CR distribution is quite flat. Thus a fiat gamma- 
ray surface brightness in the Fermi bubbles in the jet 
scenario with viscosity is not sensitive to the exact time 
of observation. 

Figure [9] also clearly shows that the CR bubble expands 
as it rises. The maximum expansion speed happens at 
the highest latitude, where it is around 2000 km/s {r^ 10 
degree/Myr at a distance of 8.5/cos(50°) ~ 13.2 kpc 
from the solar system) at the current time t = tpermi 
(1.49 Myr for run V3-diff3). At lower latitudes, the bub- 
ble surface expands at slightly lower velocities (~ 1500 
- 2000 km/s), which incline toward the original jet di- 



rection instead of perpendicular to the bubble surface 
(see velocities superposed in the top panels). Another 
interesting result shown in Figure [9] is the morphological 
evolution of the Fermi bubbles seen in Galactic coordi- 
nates. At early times, the bubbles are elongated in the 
jet (vertical) direction, but the expanding bubbles be- 
come more and more spherical and even elongated in the 
horizontal direction at later times. This is mainly be- 
cause the top and bottom of the bubbles rise in the same 
direction, while the sides of the bubbles expand in oppo- 
site directions (see top panels in Fig. [9]), and hence have 
larger relative velocities. The projection effect, i.e., the 
distance from the solar system to the bubble top is much 
larger than that to the bubble sides, also contributes to 
this morphological evolution in the Galactic coordinate 
system. 

4. SUMMARY AND IMPLICATION 

The Fermi Gamma-ray Space Telescope recently de- 
tected two extended kpc-sized gamma ray bubbles in the 
inner Galaxy. In a companion paper (Paper I) we showed 
that a recent AGN jet event from the GC could produce 
these two bubbles with roughly the observed age, loca- 
tion, size and shape. However, detailed comparisons of 
our model with observations also reveal two discrepan- 
cies: surface irregularities and non-uniform gamma ray 
limb darkening. These shortcomings motivate us to fur- 
ther explore the jet scenario for the Fermi bubbles in the 
current paper, investigating the potential roles of shear 
viscosity and CR diffusion on the bubble evolution. 

In this paper, we performed a series of hydrodynami- 
cal jet simulations incorporating shear viscosity and CR 
physics (CR dynamics, advection and diffusion). In our 
model, the opposing jets contain two components - ther- 
mal gas and CRs - and naturally form two symmetric 
CR- filled bubbles in the inner Galactic halo. As the level 
of viscosity increases, it suppresses shear motions in the 
bubble. When the viscosity coefficient is larger than a 
lower limit (/ivisc > Mmin 3 g cm"-"^ where /imin is 
even less than 1% of the Spitzer viscosity of the shock- 
heated keV-temperature gas), shear viscosity effectively 
suppresses the growth of KH instabilities at the bubble 
surface, resulting in smooth bubble edges as observed. 

The suppression of shear motions due to viscosity re- 
duces the level of CR advection, significantly affecting 
the CR distribution inside the Fermi bubbles. In non- 
viscous runs, the CR-carrying jet materials form an ax- 
isymmetric backflow along the bubble surface, which 
reaches the bubble bottom and then rises along the jet 
direction again, forming circulating motions in the bub- 
ble interior. The circulating motions mix CRs inside 
the bubble and form a CR- filled Fermi bubble. Such a 
roughly-uniform CR distribution would produce a limb- 
darkened gamma ray intensity distribution in the line- 
of-sight projected Galactic coordinate system, which is 
inconsistent with observations. However, in viscous runs 
where viscosity suppresses KH instabilities, the jet back- 
flow slows as it flows along the surface of the expanding 
bubbles. Thus the CR-carrying backflow mainly stays at 
the bubble surface, resulting in a CR distribution that 
increases toward the bubble edges, which is a key require- 
ment to produce the observed flat gamma-ray surface 
brightness distribution. Such an edge-favored CR distri- 
bution is hard to achieve by other physical mechanisms, 
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Fig. 9. — Top: Central slices (16 x 15 kpc) of log(ec) in run V3-diff3 at t = 0.7, 1, and 2 Myr. Horizontal and vertical axes refer to R and 
z respectively, labeled in kpc. Arrows superposed show thermal gas velocity. Bottom: The corresponding distributions of the line-of-sight 
projected CR energy density in run V3-diff3 at these three times in Galactic coordinates. 



suggesting that viscosity may indeed play a significant 
role in jet evolution. 

The observed Fermi bubbles have very sharp gamma- 
ray edges, indicating that CR diffusion across bubble 
edges is suppressed significantly below the CR diffusion 
rate estimated in the solar vicinity. However, if CR dif- 
fusion is also suppressed in the bubble interior, viscous 
runs produce limb-brightened gamma-ray bubbles in the 
projected Galactic coordinate system, which are also in- 
consistent with observations. Interestingly, if CR diffu- 
sivity in the bubble interior is close to that estimated 
in the solar vicinity, CR diffusion transports CRs from 
the bubble edges to the interior, resulting in a roughly 
fiat projected CR distribution. Thus, by including CR 
diffusion within, but not across, the bubble edges, our 
viscous runs produce CR-filled bubbles very similar to 
those observed. 

Given that we already invoke magnetic draping-which 
is inevitable and has been observed in numerous sys- 
tems ranging from comets to the Sun's coronal mass 
ejections — to suppress CR diffusion, it might seem su- 
perfluous to invoke viscosity as a stabilizing mechanism, 
since magnetic t ension in the drape may also stabilize the 
KH instabilitv (fLyutikov J006; Ruszkowski et al. 2007; 
iDursi fc Pfrommerll2008l ). As suming that the results o f 
subsonic draping described in iDursi fc Pfrommerl ()2008l ) 
apply here, the field strength in the drape is indepen- 
dent of the ambient galactic field strength and will build 
up to be in rough equipartition with ram pressure (for 
the parameters we have chosen, it is ~ 10 /iC); this 
means that the Alfven speed is of order the shear ve- 
locity, which is the requirem ent for the KH instability 
to be quenched (jPursil l2QQ7l ). We have four comments 
on this: (1) a further requirement for stabilization by 
magnetic tension is that the swept up ambient field has 



a coherence length Ab which is larger than or of or- 
der the bubble size R; otherwise, the bubble will still 
be sh redded apart, as demonstrated in MHD simula - 
tions (|Ruszkowski et al.ll2007l : IDursi fc Pfrommerll2008f ). 
While Ab ^ is plausible for bubbles in galaxy clus- 
ters, it is as yet unclear whether the halo of our Galaxy 
has a field coherent on ~ 10 kpc scales, without sig- 
nificant small scale irregularities. If other constraints 
on magnetic coherence ca n be placed, such as polarized 
emission from the drape (jPfrommer fc Jonathan Dursil 
I2OIOI ). this would indirectly probe the halo magnetic field 
(while there has been no detection of polarization of the 
haze/bubbles by WMAP, this niay be due to the sig- 
nificant noise in the data: lDoblerll2012l ). (2) Viscosity is 
also required to stabilize internal flows within the bubble, 
which also induce KH instabilities — indeed, our studies 
suggest that this source of shear is the dominant com- 
ponent. The influence of the swept up magnetic sheath 
is less clear in this case. (3) As discussed earlier in this 
section, viscosity effectively suppresses the backward mo- 
tions of the jet backflow and the circulating motions in 
the bubble interior, naturally leading to an edge-favored 
CR distribution, which is a key ingredient to produce the 
observed flat gamma ray intensity distribution. (4) The 
dynamics of draping in a supersonic flow where a strong 
shock forms - as simulated here - could be substantially 
different from the subsonic case simulated by previous 
authors. 

We also studied the evolution of the Fermi bubbles 
in our simulations. The top and side of the bubbles 
are expanding explosively, currently at speeds of around 
1500 — 2000 km/s, while the bottom is rising along the 
jet direction. In the Galactic coordinate system, the bub- 
bles are elongated along the jet direction at early times 
(and the current time), and then become more spherical 
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and even elongated in the perpendicular direction at late 
times. The flatness of the projected gamma ray intensity 
distribution may last for a quite long time (> 1 Myr). 

This work strengthens the previous result in Paper I 
that the Fermi bubbles may have been recently produced 
by two opposing CR-carrying jets. The true level of vis- 
cosity in a magnetized low-collisional hot plasma is un- 
clear. In particular, nearly-parallel magnetic field lines 
at the bubble surface may significantly suppress momen- 
tum transport across bubble edges. The level of viscosity 
in collisionless plasma in the bubble interior is even more 
uncertain. However, our calculations indicate that even 
a significantly-suppressed low level of viscosity (~ 0.1% 
- 1% of the Spitzer viscosity in the shocked surrounding 
gas) can play a significant role during the jet evolution, 
producing smooth bubble edges and an edge-favored CR 
distribution, which are difficult to accomplish by other 
mechanisms. If momentum transport across the bubble 
surface is fully suppressed by magnetic fields, we specu- 
late that viscosity in the bubble interior alone may be 
sufficient to suppress the backward motion of the jet 
backflow. 



This paper suggests the potentially important roles of 
viscosity and CR diffusion in the Fermi bubble event; our 
simulations assume these to be isotropic, as is perhaps 
appropriate for highly tangled fields. To advance these 
suggestions further requires MHD simulations which are 
able to take into account the anisotropic nature of viscos- 
ity and CR diffusion, as well as other effects such as mag- 
netic draping. We note that our simulations nonetheless 
lay the groundwork for establishing that physically rea- 
sonable values of viscosity and CR diffusion-which are 
very uncertain — can help explain the main features of 
the Fermi bubbles. 
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APPENDIX 

IMPLEMENTATION OF COMPRESSIBLE VISCOSITY 

In this Appendix, we explicitly present the compressible viscous terms in CR- hydro equations [2] and [3] and briefly 
show how they are implemented in our code. We only consider shear viscosity and neglect bulk viscosity. In Cartesian 
coordinate systems, the viscous stress tensor 11 can be written as 

However, in the cylindrical coordinate system (i?, 6, z), the expression for 11 is much more complicated. In this paper, 
we assume axisymmetry: 

^=0 and ^^ = 0, (A2) 
which significantly simplifies the formulae for the components of 11. 
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We first calculate the tensor of the velocity gradient Vv in physical units: 

Vv 



'dvn/dR 
vr/R 
dvJdR 



dvYi/dz^ 


dvz/dz 



(A3) 



and define the tensor A: A = [Vv + (Vv)*^]/2, where (Vv)*^ is the transpose of Vv. Then the viscous stress tensor 
can be written as (in physical units): 

n = 2/ivisc[A-i(V-v)I], (A4) 
where I is the identity matrix Jij = (^ij, and the divergence of gas velocity is 

After calculating all non-zero components of 11 for each simulation zone, we finally calculate the viscous term V • 11 
in the gas momentum equation [2] and 11 : Vv in the gas energy equation [3l 



V n 



^ + ^(nRR-n,,) + ^,o, ^ + — + ^ 



n : Vv = n : A = ^ Eij 



Since viscosity is explicitly implemented in our code, the time-step is constrained to ensure numerical stability: 

?2 



dt < dt 



vise — ^VISC 



Cvisc min 



dUt dzf 



(A6) 
(A7) 

(A8) 



where pij is the gas density in the grid cell zj), the constant Cvisc is chosen to be Cvisc = 0-3, and the minimization 
occurs over all grid zones. 



